# -*- coding: utf-8 -*-
"""
Created on Thu Apr 23 08:43:36 2015

@author: ba0hx
"""

import numpy as np
import scipy.stats as ss


def case(n = 10, mu = 3, sigma = np.sqrt(5), p = 0.025, rep = 100):
    m = np.zeros((rep, 4))


    for i in range(rep):
        norm = np.random.normal(loc = mu, scale = sigma, size = n)
        xbar = np.mean(norm)
        low = xbar - ss.norm.ppf(q = 1 - p) * (sigma / np.sqrt(n))
        up = xbar + ss.norm.ppf(q = 1 - p) * (sigma / np.sqrt(n))

        if (mu > low) & (mu < up):
            rem = 1
        else:
            rem = 0

        m[i, :] = [xbar, low, up, rem]

    inside = np.sum(m[:, 3])
    per = inside / rep
    desc = "There are " + str(inside) + " confidence intervals that contain "
    "the true mean (" + str(mu) + "), that is " + str(per) + " percent of the total CIs"


    return {"Matrix": m, "Decision": desc}

def case2(n = 10, mu = 3, sigma = np.sqrt(5), p = 0.025, rep = 100):
    scaled_crit = ss.norm.ppf(q = 1 - p) * (sigma / np.sqrt(n))
    norm = np.random.normal(loc = mu, scale = sigma, size = (rep, n))

    xbar = norm.mean(1)
    low = xbar - scaled_crit
    up = xbar + scaled_crit

    rem = (mu > low) & (mu < up)
    m = np.c_[xbar, low, up, rem]

    inside = np.sum(m[:, 3])
    per = inside / rep
    desc = "There are " + str(inside) + " confidence intervals that contain "
    "the true mean (" + str(mu) + "), that is " + str(per) + " percent of the total CIs"
    return {"Matrix": m, "Decision": desc}

print case2()["Decision"]
print case2()["Matrix"]